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QCD is expected to have a rich phase structure. It is empirically known to be difficult to 
access low temperature and nonzero chemical potential fi regions in lattice QCD simulations. 
We address this issue in a lattice QCD with the use of a dimensional reduction formula of 
the fermion determinant. 

We investigate spectral properties of a reduced matrix of the reduction formula. Lattice 
simulations with different lattice sizes show that the eigenvalues of the reduced matrix follow 
a scaling law for the temporal size Nt- The properties of the fermion determinant are 
examined using the reduction formula. We find that as a consequence of the Nt scaling 
law, the fermion determinant becomes insensitive to /i as T decreases, and /i-independent at 
T = for /i < 

The Nt scaling law provides two types of the low temperature limit of the fermion 
determinant: (i) for low density and (ii) for high-density. The fermion determinant becomes 
real and the theory is free from the sign problem in both cases. In case of (ii), QCD 
approaches to a theory, where quarks interact only in spatial directions, and gluons interact 
via the ordinary Yang-Mills action. The partition function becomes exactly Z-^ invariant 
even in the presence of dynamical quarks because of the absence of the temporal interaction 
of quarks. 

The reduction formula is also applied to the canonical formalism and Lee- Yang zero 
theorem. We find characteristic temperature dependences of the canonical distribution and 
of Lee- Yang zero trajectory. Using an assumption on the canonical partition function, we 
discuss physical meaning of those temperature dependences and show that the change of the 
canonical distribution and Lee- Yang zero trajectory are related to the existence/absence of 
^-induced phase transitions. 
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§1. Introduction 

QCD, first principle of strong interaction, has two important features; confine- 
ment and chiral symmetry breaking. They change its nature depending on circum- 
stances, e.g., temperature (T) and quark chemical potential (/i), which leads to rich 
structure in the QCD phase diagram.^®'® Hadrons and nuclei are formed at or- 
dinary temperatures and chemical potentials. The quark-gluon plasma (QGP) is 
formed at high T, which is expected to be created in the early universe and also 
in experiments, RHIC, LHC, etc. Several possibilities have been proposed for high 
density states of matter, which are realized in the core of neutron stars. 

Finite temperature properties of QCD have been uncovered by lattice QCD sim- 
ulations ® which is a computational approach implemented with Monte 
Carlo methods. For instance, recent simulations on finer lattices show that the de- 
confinement transition is crossover,® and occurs at Tpc = 150 — 170 MeV depending 
on observables.^*'^''"^' 

In contrast, the properties of QCD at nonzero /i have been difficult to study 
because of the notorious sign problem.^® The sign problem spoils the importance 
sampling at fj,{^ 0), and makes it unfeasible to generate gauge ensemble. Several ap- 
proaches have been developed to avoid the sign problem and study nonzero- systems 
in lattice QCD simulations, see e.g., Refs.'m''1^''I^'n3I'[S} xhe consistency between 
different approaches are shown for small n/T. For instance, it was 

showrPJ-EJ that 

several approaches of finite density lattice studies are consistent up to ///T ~ 1, for 
staggered fermions. For Wilson fermions, we showed the consistency of the Tay- 
lor expansion, multi-parameter reweighting (MPR) and imaginary chemical poten- 
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tial methods.'^ For status of lattice studies of the QCD phase diagram, see e.g, 
RefsP'l^ Extensive studies have been made for the finite density properties of 
QCD in particular for the location of the QCD critical end point (CEP). 

The QCD phase diagram contains rich physics also at low temperaturesP''^''^ 
One expectation is that towards large density region, QCD changes its state from 
hadron gas, nuclear liquid and color superconducting state. In addition, the discovery 
of a pulser with twice solar masd^ calls the reliable equation of state based on QCD 
at low temperature and high density. The study of the low temperature and finite 
density states would be an important challenge for lattice QCD simulations. 

There are expectations of the existence of sign free regions at low temperatures. 
The finite density lattice QCD at low temperatures had been extensively investigated 
in the Glasgow method, see e.g. Ref.,'23'[22li,[23li,[24li -^^j-^gj-g -^g^g ghown the onset of 
the baryon density is given by m,r/2. Recently, the phase at low temperatures was 
investigated in a lattice simulation .^29 It is empirically expected that the fermion 
determinant is independent of /x at T = up to ;U < Mjv/3. However, the inclusion 
of chemical potentials can generally change the eigen spectrum of the Dirac operator 
and hadron spectrum. The //-independence at T = was, therefore, considered as a 
small puzzle and called "Silver Blaze" problem™-' Cohen showed the //-independence 
at T = for isospin chemical potential fij case and Adam^23 for quark chemical 
potential case. The orbifold equivalence for baryon chemical potential and isospin 
chemical potential also suggestecP^'l^ that the phase of the fermion determinant is 
small up to m,r/2. The same result was also obtained by Splittorff and Verbaarschot 
by using the chiral perturbation theory (ChPT). On the other hand, the mechanism 
to extend the Silver Blaze region from fj, = to // = Mn/3 is not yet understood. 

A sign-free region was also suggested for high density. Hong derived a high density 
effective theory (HDET) of low energy modes in dense QCD // ^ ^qcd-'^ The 
positivity of the fermion action in HDET was later discussed by Hong and Hsu.l^D'ISS 
HDET was also studied in e.g., Refs.lSf'® 

In this paper, we study the property of QCD at low temperature and finite 
density, using lattice QCD. Particularly, the properties of the fermion determinant 
at low temperature and nonzero /x will be clarified. The key idea to tackle this 
issue is a reduction formula for the fermion determinant. The formula is obtained by 
performing the temporal part of the determinant of a fermion matrix. The technique 
is analogous to the transfer matrix method in condensed matter physics, and provides 
several advantages in lattice simulations of QCD at finite density. The formula 
produces a matrix with a reduced dimension, which is called the reduced matrix. 
We will see that the reduced matrix plays an important role in the properties of the 
fermion determinant at low temperature and nonzero /x. 

Some issues are discussed in the present paper. In § [21 we focus on the reduced 
matrix, and discuss its derivation, interpretation and spectral property, where the 
connection with important dynamics of QCD is considered. In particular, we will find 
a clear indication from lattice QCD simulations that the eigenvalues of the reduced 
matrix follow a scaling law of the temporal size Nf. This section is partly a review 
of the reduction formula. 

In § [21 the property of the fermion determinant is investigated in detail. We 
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show that the quark determinant is insensitive to /x at low temperatures for smah ^u. 
Using the relation between the eigenvalues of the reduced matrix and pion mass, we 
will show that the /i-independence continues up to ^ = m7r/2. We will see that the 
Nt scaling law leads to two expression of the low temperature limit of the fermion 
determinant both for low density and for high density. The low density expression 
corresponds to the /i-independence at small n. Hence the //-independence at T = 
for fj. < niT^/l is a consequence of the A^j-scaling law of the eigenvalues of the reduced 
matrix. The fermion determinant becomes real and the theory is free from the sign 
problem both cases. 

In the case of high density and low temperature limit, QCD approaches to a 
theory, where quarks interact only in spatial direction, and the gluon action is given 
by the ordinary Yang-Mills action. The corresponding partition function is exactly 
Z3 invariant even in the presence of dynamical quarks because of the absence of 
the temporal interaction of quarks. Property, possible application and numerical 
feasibility are discussed. 

In § U the reduction formula is applied to the canonical formalism and Lee- Yang 
zero theorem. These methods are often used in analysis to find CEP and first order 
phase transition line at low temperatures. The sign problem causes difficulties in 
these methods .E3:[ISli For instance, Ejiri pointed out that the sign problem causes 
a fictitious signal for the finite size scaling analysis of the Lee- Yang zero closest to 
the positive real axis, which makes it difficult to distinguish crossover and first order 
transition. 

We consider the fugacity expansion of the grand partition function, where ZccitJ^) 
is expressed via the canonical partition functions Zn with quark number n. We show 
that both the n-dependence of and the Lee- Yang zero trajectory show character- 
istic change from high to low temperatures in a correlated manner. These behavior 
can be used to qualitatively distinguish the crossover and first order phase transition. 

§2. Reduction Formula 

In this section, we study the reduction formula of the fermion determinant. A 
basic idea of the reduction formula is to analytically carry out the temporal part 
of the fermion determinant, which reduces the rank of the determinant. Since the 
chemical potential /x couples only to the temporal link variables, the reduction for- 
mula rearranges the fermion determinant regarding /i, which enables us to separate 
out /i from link variables. These characters offer some advantages in finite density 
simulations. 

The reduction formula is not only a technical tool but also has physical interpre- 
tation. The formula produces a matrix with a reduced dimension, which we refer to 
as the reduced matrix. Its eigenvalues characterize the //-dependence of the fermion 
determinant. The reduced matrix physically corresponds to a transfer matrix or a 
generalized Polyakov line, and its eigenvalues are related to a free energy of dynam- 
ical quarks. 

In this section, the spectral properties of the reduced matrix are investigated in 
detail. The formulation is given in § 12.11 The interpretation and overview of the 
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reduced matrix are presented in § 12.21 Simulation setup is explained in § 12.31 Some 
important and fundamental properties of the reduced matrix are shown in § 12. 4[ In 
§ \2.b\ we will find a clear indication of the A^^-scaling law of the eigenvalues of the 
reduced matrix. In § 12.61 we discuss the relation between the eigenvalues of the 
reduced matrix and pion mass. 

2.1. Formulation 

The reduction formula was derived by GibbspS for staggered fermions. Alterna- 
tive derivation was found by Hasenfratz and Toussaint.'^ Later it has been applied 
to various studies, e.g. Glasgow methodpD'E2li>[23]i,[24ji,[38li nmlti-parameter reweight- 
ingpSJ canonical formalism.HS'CSli'EJ'liSl' The derivation for Wilson fermions needs a 
delicate treatment because the Wilson terms (r ±74) are singular for r = 1. Adams 
derived the formula using zeta-regularization.SS Borici derived the formula by intro- 
ducing a permutation matrix.SS Borici's method was further studied by Alexandru 
and Wengei^Si and two of present authors.l^SJ The formula for the Wilson fermion 
was later applied to studies of CEpf^S and thermodynamical quantities.'^ The re- 
duction formula in continuum theory was found by Adams using zeta-regularization, 
which was used to show the /x-independence of the fermion determinant at T = 

Throughout the present paper, we consider the Wilson fermions. For the deriva- 
tion for the staggered fermions, see e.g. Ref.^ The grand partition function of 
A^j-flavor QCD at a temperature T and quark chemical potential /i is given by 

Zgc{i^,T) = j^U [detZ\(/x)]^/e-^«. (2-1) 

In our simulations, a renormalization group(RG)-improved action is used for Sq- 
The clover-improved Wilson fermion is employed for the fermion action 

3 

1=1 

— i^Csw^x,xi ^ (y ii.vFii.v, (2-2) 

where k and r are the hopping parameter and Wilson parameter, respectively. The 
quark chemical potential /x couples to the fourth current V'74V' in Euclidean path- 
integral, then a temporal hop accompanies a factor e^^'^W Csw 

is a coefficient of 

the clover term. Note that the clover term is diagonal in the coordinate space, hence 
it does not cause any problem in the derivation of the reduction formula. 
The fermion determinant satisfies a 75-hermiticity relation 

(det Z\(/i))* = det Z\(-/x*). (2-3) 

It ensures that det Z\(0) G M. In the presence of nonzero real /i, the fermion determi- 
nant is in general complex, which causes the sign problem. The det A is real if /i is 
pure imaginary. This property has been used in the study of the QCD phase diagram 
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in lattice QCD simulations^ ^'^'^^^^^'^'^^'^^'^'^ and also in 
phenomenological studies, see e.g. Refs.'^''^ 

For the preparation of the reduction formula, we define the spatial part of the 
Wilson fermion matrix as B, 

3 
i=l 

— nCswSx,x' E cr^i^F^t,, (2-4a) 

IJ,<U 

and two block-matrices 

= c-B^^^f'^ix, y, ti) r'L" - 2c+k r^5'^H[x - y), (2-4b) 
P, = (3-'''^-{x,y,ti) 

= c+B'^^'^'^x, y, ti) rl^Ul\y, U) - 2c^k r'^J' 5{x - y)Uf{y, U)- (2-4c) 

Here r± = (rib 74)72, which are projection operators in case of r = 1. This property 
is used in the derivation of the formula. c± are arbitrary scalar except for zero, which 
is introduced in so called permutation matrix P = (c_r_ + c+r+V^e'^"). Oi describes 
a spatial hopping on a fixed time plane ti, while /3j describes a spatial hopping at ti 
and a temporal hopping to the next time slice. They are independent of /x. 

A temporal matrix representation of the Wilson fermion matrix contains block- 
elements proportional to r±, which are singular. This fictitious singularity is avoided 
by calculating det{PA)^^ Then the determinant is obtained from det(PZ\)/ detP. 
Using the block matrices, the reduction formula is given bjB^ 



detZ\(;u) = {c+c^)-^/'^C^'^''^^Codet + Q) , (2-5a) 

with 

Q = (a-i/3i)...(a^i/3^J, (2-5b) 

Co = (fldetiai)] , (2-5c) 



where ^ = exp(-/i/T) is the fugacity. N = iNcN^Nt and A^red = N/Nt are the 
dimension of A and Q, respectively. The rank of and Q is reduced by 1/Nt 
compared to A^. Instead, Q is a dense matrix. Q and Cq are independent of ^. 
Hence, the reduction formula separates the chemical potential from the gauge parts. 

To obtain detZ\, we need to evaluate det{Q + S,)- Here we calculate the eigen- 
values A for IQ — A/| = 0. Once we obtain A, the quark determinant is the analytic 
function of fi. Then, values of det A{fi) are obtained for arbitrary values of fi, which 
is an advantage, although the eigen problem requires large numerical cost. Alterna- 
tive methods such as LU decomposition of Q + ^ are available instead of solving the 



'Towards extremely dense matter on the lattice 



7 



eigenvalue problem for Q. In this case, we need to perform the LU decomposition 
for each value of /i. 

Once having the eigenvalues of Q, we obtain 

det Aip) = CoC^^'^''^^ n + ^) (2-6a) 

n=l 

= Co X; Cnr-^-<^/' = Co E c-^"' (2-6b) 

n=0 n=-Ar,ed/2 

where we set c± = 1 for simplicity. In the second line, we redefine the index n from 
the first expression to the second one. The fermion determinant is described in two 
forms: a product form Eq. (12 -Gap , and a summation form Eq. (|2-6bp . The second one 
is nothing but a fugacity or winding number expansion of the quark determinant, 
where fugacity coefficients c„ are polynomials of A„.^^' The consistency between 
the reduction formula for the Wilson fermions and staggered fermions can be found 
by considering the fugacity expansion form. The fugacity expansion of the fermion 
determinant is also obtained in other approaches, such as a domain decomposition 
technique.'^ 

Here, we summarize some features of the reduction formula. The formula offers 
several advantages; 

• The dimension of the determinant is reduced by 1/Nt, which reduces the nu- 
merical cost of the direct calculation of the fermion determinant. 

• det A{fj,) is the analytic function of fi, which suppresses the numerical cost to 
evaluate values of det /\(//) for many values of /i and also provides an insight 
into the //-dependence of det Z\(/u). 

• Increase of the numerical cost for low temperature(large Nt) is also suppressed. 
On the other hand, it has an applicable range; 

• A direct method is used for the eigenproblem of the reduced matrix. This 
requires large numerical cost. In particular, it is difficult to apply the reduction 
formula to large lattice size because of the limitation of a memory. 

2.2. Physical interpretation 

In the reduced matrix Q, block elements (a~^/3j) describe propagations of a 
quark from t = ti to t = tj+i. Q is the product of the block elements, and therefore 
it describes propagations from the initial t = ti to final time t = tNt (depicted in the 
right panel of Fig. [T| . Accordingly, the reduction formula is analogous to the transfer 
matrix method, where block elements is interpreted as transfer matrices [3Sli:lli)>ll3 
The matrix Q is also understood as a generalized Polyakov line. If we pick up the 
temporal links in Q, it is written as Q = ■ ■ ■ U4{ti) ■ ■ ■ U4^{t2) ■ ■ ■ (/^{tN^). If the trace 
is taken, trQ is similar to the Polyakov loop. It is known that the Polyakov loop 
describes the free energy of static quarks (P) ~ e~^^'^ with P = tr J^^j^ f74(tj). 
However, Q contains spatial hopping terms denoted as the dots, where quarks are 
not static. This suggests that the reduced matrix is related to the free energy of 
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Fig. 1. Schematic figures for Co and Q. The left panel describes spatial 
loops on a plane t — ti. Co consists of those spatial loops. The right 
panel denotes Q, which describes propagations of quarks from t = ti 

to t = tNf 



't = t 




Fig. 2. Schematic figures for the reduction procedure, where the spatial 
loops Co are separated from Q. 



dynamical quarks. 

Co consists of the spatial loops, where each loop is in a fixed time slice, which 
accompanies no temporal hops, see the left panel of Fig. [1] and Fig. [2j Thus the 
spatial quark loops Co are separated from temporal hoppings and do not contribute 
to the /^-dependence of the fermion determinant. 

Here, we summarize properties of the reduction formula. From 75-hermiticity, 
it follows ; 

• Co is real. 

• Two eigenvalues form a pair 

Xn,l/Xn- (2-7a) 

This is because of a symmetry of Q, see e.g., Ref.^ We give a simple proof of 
this relation in Appendix. lA.li 

• The coefficients of the positive and negative winding terms satisfy 

C-n = <• (2-7b) 

• The product of all the eigenvalues is unity (followed from c_„ = c*); 

det Q = JJ A„ = A1A2 • • • A^^^, = 1. (2-7c) 

n=l 
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• For the case c± = 1, it is straightforward to show 

deta, = det/3i. (2-7d) 

• From Eq. ()2-7ap and Eq. ()2-7cp . we can separate all eigenvalues into two sets 
({'^fcll-^fcl > 1} and {A/fc||Afc| < 1} ). Then, the product of the normalized 
eigenvalues in each set is ±1, 

For later convenience, we introduce a notation for |A| > 1 and Ag for |A| < 1. 
These properties are satisfied for configuration by configuration, independent of 
temperature. 

Although it is non-trivial to clarify the correspondence between the reduced 
matrix Q and the original matrix Z\, properties of QCD can be seen in the spectral 
property of Q; e.g., 

Confinement: As we have mentioned, the matrix Q is related to the Polyakov 
line. A correlation between the Polyakov loop and the eigenvalues of the reduced 
matrix was found in Ref.l^ We can find the confinement property is realized in 
the angular distribution of A. 

Chiral symmetry breaking : It is known that the distribution of the eigenvalues 
form a gap near |A| = 1. Gibbs pointed oulPS that the gap of the eigenvalue 
distribution is related to the pion mass. Hence, the behavior of the eigenvalues 
near the unit circle is related to the chiral symmetry breaking. 

Hadron Spectroscopy: Since the reduced matrix is related to the free energy of 
dynamical quarks, its eigenvalues are related to hadron spectrum. Fodor, Szabo 
and Toth proposed to extract hadron masses from the eigenvalues of the reduced 
matrix based on thermodynamical approach.'^ 

Low temperature behavior : Because Q is a product of Nt block-matrices, then 
it is expected^' that the eigenvalues of Q follow a scaling law for the temporal 
lattice size Nt- The scaling law is useful to study the low T behavior of the 
fermion determinant. For instance, we will see that the scaling law explains the 
/^-independence of the fermion determinant at T = for fi < 

2.3. Simulation setup 

Before proceed to numerical simulations, simulation details are given here. Sim- 
ulations were performed on four lattice setups; (N^ x Nt, mps/my) = (i) (8^ x 4, 0.8) 
(ii) (8-^,0.8), (iii) (10^ x4,0.8) and (iv) (8^ x4,0.6). The wide range of temperatures 
were covered by using (i) and (ii). The setup (iii) and (iv) were used to investigate 
the finite size effect and quark mass dependence. The detail of these simulations is 
as follows. 

(i) We investigated 29 values of /? in the interval 1.5 < /? < 2.4, which covers 
the temperature range T/Tpc = 0.76 — 3. The data were used in our previ- 
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ous studies on imaginary chemical potential approach^^' and thermodynamical 
quantities!^ 

(ii) This simulation was performed to examine lower temperature. We investigated 
16 values of /3 in the interval 1.79 < (3 < 1.94 corresponding to the temperature 
range T/Tpc = 0.460 - 0.587. 

(iii) This setup was performed to study the finite size effects. We investigated 29 
values of /3 in the interval 1.79 < (3 < 1.94, which covers the temperature range 
T/Tpc = 0.93 - 1.2. 

(iv) This simulation was performed to investigate the quark mass dependence. We 
investigated four values of /3 = 1.6,1.7,1.95 and 2.0, which correspond to 
T/Tpc = 0.86,0.94, 1.48 and 1.67, respectively. 

For all the four setup, the value of the hopping parameter k was determined for 
each /3 by following lines of constant physics (LCP) in Ref.-^ The clover coefficient 
Csw was determined by using a result obtained in the one-loop perturbation theory 
: Csw = (1 - 0.8412/3"^ We also used the data in Ref.E^ to obtain the values 
of the temperature from (3. 

Gauge configurations were generated at fi = with the hybrid Monte Carlo 
simulations. The size of the molecular dynamics step was St = 0.02 for mps/w-v = 
0.8 and 5t = 0.01 for mps/my = 0.6 with unit length. The acceptance ratio was 
more than 90 % for Ng = 8 and 80 % for A^^ = 10. HMC simulations were carried 
out for 11, 000 trajectories for each parameter set for mps/my = 0.8, and 4, 000 
trajectories mps/w-v = 0.6. It should be noted that the statistics are small for 
mps/mv = 0.6. 

The measurements were performed for each 20 HMC trajectories after thermal- 
ization. The eigenvalue calculations were performed for 400 configurations for heavy 
quark case, and 50 configurations for light quark case. Computational details of the 
eigenvalue calculation were explained in Ref.l^ 

2.4. Pair nature, Z3 symmetry and temperature dependence 

Let us study the spectral properties of the eigenvalues of the reduced matrix. 

We look at typical behaviors of the eigenvalue distribution on the complex A 
plane in Fig. [3l where results are shown for three temperatures. Each result is 
obtained from one configuration on the 8'^ x 4 lattice. For each temperature, the 
eigenvalues are shown in two different scales; the top panels show the distribution of 
whole the eigenvalues in the complex plane, and bottom panels enlarge the region 
near the origin. These panels show two features of the reduced matrix, the pair 
nature and Zs-like property. 

First, we consider the angular (phase) distribution of the eigenvalues of the 
reduced matrix. The histogram of the angular distribution is shown in the right panel 
of Fig. m We find that the histogram of the angular distribution is well fitted with the 
three Gaussian function p{9) = J2i=i 23^* exp(— (6* — Cj)^/(2c^)), where = arg(A). 
Three peaks are located at 6 = 0, ±2^/3 at T < Tpc, and two of them (9 = ±2^/3) 
shift towards ^ = as T increases. Although Z3 is explicitly broken in the presence 
of the quarks, the eigenvalue distribution is correlated to the Polyakov loop, as we 
have mentioned in the previous subsection. Indeed, Alexandra and Wenger found 
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Fig. 3. Eigenvalue distribution on the complex A plane for three temperatures on the 8'^ x 4 lattice 
with mps/mv = 0.8. The top panels show the distribution of whole the eigenvalues, while the 
bottom panels magnify the region near the origin of the top panels. 
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Fig. 4. The spectral density of the larger half of the eigenvalues An for three temperatures. The 
result for the smaller half of the eigenvalues is obtained from the pair nature. Left : absolute 
value, right : phase of A„. In the right panel, thick lines are obtained from the three-Gaussian 



function : p{6) ; 
arg(A) = ±27r/3. 



2 3 ^^P(~(^ — Ci) /(2Ci)), {9 = arg(A)). The vertical lines denote 



a correlation between the argument of the Polyakov loop and the location of the 
blank region of the eigenvalue distribution at high The confinement property 

of QCD may be seen in the angular distribution of the eigenvalues of the reduced 
matrix. The left panel of Fig. U] shows the distribution of the magnitude of the large 
eigenvalues A^. The distribution of \Xl\ broadens at high temperatures. 

Next, we consider a gap of the eigenvalues. According to the pair nature, the 
eigenvalues are split into two regions. The half of the eigenvalues are distributed in 
a region |A| < 1, and the other half in a region |A| > 1. There is the gap between 
the two regions where no eigenvalue exists, see the bottom panels of Fig. [21 As we 
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will see later, the gap is related to the pion mass. 

Let us consider a physical meaning of the two types of the eigenvalues: larger 
half (|A| > 1) and smaller half (|A| < 1). The pair nature of the eigenvalues is a 
consequence of 75-hermiticity, and therefore of the symmetry between the quark and 
anti-quark. As we have mentioned, tr Q is related to a generalized Polyakov loop 
and to free energies of an quark and anti-quark. Qualitatively, the relation can be 
written as A ~ e~^^'^ with a free energy F. Obviously [A| < 1 if F > 0, while |A| > 1 
if -F < 0. Hence, the smaller half of the eigenvalues correspond to quarks and the 
larger half correspond to anti-quarks. 



Q 




Q 



800 




-3-2-10123 
Phase (k) 



Fig. 5. Finite size effect on the spectral densities of tfie larger half of the eigenvalues. Left : absolute 
value, right: phase. In each panel, the blue-solid and red-dotted lines denote the results for 8^ x 4 
and 10=^ X 4. Note Vs{Ns = 10)/K(iV, = 8) ~ 2. 

The volume dependence of the eigenvalues is shown in Fig. [SI It turns out that 
the spectral density of the eigenvalues is almost insensitive to the spatial lattice 
size Ng. The effect of increasing Ng is small even for the gap and tail part of the 
distribution. This insensitivity may be a consequence that long-range correlations are 
suppressed due to the heavy quark mass used in this work. The volume dependence 
should be investigated in future simulations both with small quark masses and large 
lattice sizes. 

Although the histogram of eigenvalues of Q is insensitive to Ns, the fermion 
determinant depends on A^^. The fermion determinant can be written as det ~ 



^^P(X]r),=i + 0)1 spectral representation, 



det A{fi) = Co exp (^2NcVsfi/T + ANcVs j dXp{\) ln(A + , 



(2-8) 



where Vg = N^. p{\) is the spectral density on the complex A plane. The term 
ln(A + ^) is complex, and generates the phase of detZ\(^). As we have shown in 
Fig. O p{\) is insensitive to Ns, and therefore the A-integral is also insensitive to Ng- 
The imaginary part of the A-integral, / dXp{\)\n.{\ + ^), is also insensitive to A^^- 
This means the phase of det A{p) is proportional to Vg, which shows the well-known 
faclP that the sign problem becomes severe for large lattice volume. 

Note that Nt dependence is included in A itself, but does not appear in the 
overall factor of the exponent. Hence, the increase of Nt does not increase the phase 



'Towards extremely dense matter on the lattice 



13 



of the determinant. This suggests that the sign problem is not necessarily severe 
at low temperatures. We will see later, the increase of makes the sign problem 
milder. 



0.03 



T/Tpc=0.93-0.94 




1000 



0.03 



T/Tpc=1.67-1.69 



m /mv= 
mp3/mu=0.8 




1000 



Fig. 6. The quark mass dependence of the spectral density. Here we consider the absolute values 
of the larger eigenvalues |A| > 1. The parameter is /? = 1.7 for rrips/mv = 0.6 and /3 = 1.8 for 
TTips/mv ~ 0.8. The decrease of the quark mass narrows the width of the distribution. 
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Fig. 7. The quark mass dependence of the scatter plot of the eigenvalues. The eigenvalues near the 
unit circle, which is related to the pion mass, is shown. Eigenvalues approach to the unit circle 
for small quark mass. 



The quark mass dependence of the eigenvalues of the reduced matrix is shown 
in Figs. Eland [71 The eigenvalues A depend on the quark mass for both the tail part 
and gap part. It is shown in Fig. [6l that the decrease of the quark mass narrows 
the eigenvalue distribution. Figure [3 shows the quark mass dependence of eigen- 
values near the unit circle. Qualitatively, eigenvalues approach to the unit circle 
as the quark mass becomes smaller. Gibbs pointed oulP^ for staggered fermions 
that eigenvalues outside the unit circle move away from the unit circle as the quark 
mass increases. The results in Figs. El and [3 are consistent with the result in Ref.'® 
However, the quark mass dependence is considered without the ensemble average. 
We will consider the quark mass dependence of the average of the gap later. 
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2.5. Low-T behavior and Nt-scaling law 



I n;:' P=1 -86 : 

k l\ 

-1.5 -1 -0.5 0.5 1 1.5 2 

In \X\/N^ 

Fig. 8. Left : The distribution of the large eigenvalues for /? — 1.8 on the 8* lattice. See also 
the result on 8^ x 4 in Fig. [3] Right : Histogram of the eigenvalue distribution with scaled by 
Nt- The blue and red lines denote the results for Nt = A and 8, respectively. The agreement 
between them implies the scaling law A ~ Z^' . The two peaks are a consequence of the pair 
nature A and 1/A*. 

Next, we investigate the Nt dependence and the properties at lower temperature. 
The left panel of Fig. [8] shows the eigenvalue distribution for the 8^ lattice with 
/3 = 1.8, which corresponds to T/Tpc ~ 0.5. Since T = l/{aNt), the temperature in 
Fig.[8]is almost half compared to the case of Fig[3l A major difference between Nt = 4: 
and A't = 8 is the magnitude of the eigenvalues. As the temperature decreases, the 
smaller half of the eigenvalues become smaller and larger half of the eigenvalues 
become larger. 

The right panel of Fig. [8] shows the histogram of the absolute value of A scaled 
by Nt- The results for A'^t = 4 and Nt = 8 agree well. This agreement implies that 
the eigenvalue distribution as a function of a variable In |A|/At is almost independent 
of Nt, which leads to the scaling law |A| ~ with a quantity In / = In |A|/A(. Note 
that / depends on the lattice spacing a. We have pointed out this scaling law in 
the previous study.^^ The matrix Q is given as the product of Ai = a~^/3i,{i = 
1, - ■ ■ Nt). Let us express Ai = A + 5Ai, where ^ is a matrix independent of time. 
Ai is expected to depend on time moderately in equilibrium, i.e., \Ai\ > \5Ai\. An 
additional cancellation would occur in ^iSAi, since SAi is a fluctuation from A. If 
this argument holds, Q is parameterize as Q ~ A^^. The agreement observed in 
Fig. [8] clearly indicates this scaling behavior. 

The right panel of Fig. [8] also shows the gap and pair nature of the eigenvalues. 
We focus on the gap in the next subsection. 

2.6. Gap and pion mass 

As we have mentioned, the reduced matrix describes the generalized Polyakov 
line and is related to the free energy of a dynamical quark. Combining several 
Polyakov lines, it is possible to construct states having quantum numbers of hadrons. 
Thus it is natural to expect that the eigenvalues of the reduced matrix have some- 
thing to do with the hadron spectrum. Indeed, Gibbs pointed oulPSJ a relation at 
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T = between the pion mass and the eigenvalue gap. Later, Fodor, Szabo and 
Toth related the hadron spectrum to the eigenvalues of the reduced matrix, based 
on thermodynamical consideration.ESJ In this subsection, we study the pion mass 
using an eigenvalue close to |A| = 1, which controls the gap. 
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T/Toc=0.93 
T/T^c=1-08 
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Fig. 9. Left panel: Eigenvalues near the unit circle |A| — 1. The parameters for each result are 
iP, Nt) = (1.8, 8), (1.9, 8), (1.8, 4), and (1.9, 4) for T/Tpc = 0.47, 0.54, 0.93, and 1.08, respectively. 
Each result is obtained from one configuration. Right panel : The largest eigenvalue in smaller 
half max|;^^|<l |A„|, which is the eigenvalue closest to the unity |A| = 1. The gap is defined as 
the deviation between the data and one. 

The temperature dependence of the gap is shown Fig. [9l The left panel shows 
twenty eigenvalues close to unity, where the result for each temperature is obtained 
from one configuration. It is shown that the gap shrinks as the temperature increases. 
To take into account the ensemble average, we focus on the maximum eigenvalue 
among the smaller half: max^x\<i ^- The gap is given by the difference between 
unity and this quantity. The result is shown in the right panel, where the gap is 
large at low T and decreases as T increases. In case of mps/mv = 0.8, the gap is 
saturated at high T. 

The right panel of Fig. [9] contains the results for two values of mps/my- We 
find that the average value of max^x\<i \M larger for smaller quark mass. This 
implies that the gap shrinks for smaller quark mass, which is consistent with Ref.'S^ 
However, the simulations for mps/my = 0.6 were done only for four temperatures 
and for the small number of HMC trajectories. Further simulations are needed to 
reveal the quark mass dependence of the eigenvalues of the reduced matrix. 

According to Gibbs,!^ the pion mass is given by 

I 

arriTr = max InjAfcp, (2-9) 

Nt |Afc|<l 

which is expected to be valid at T = 0.122} Qn the other hand, Fodor, Szabo and 
Toth derived a modified expression,!® 



am. = Jun^[-^^ln(\Y^X,\")). (2-10) 
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Fig. 10. Pion mass obtained from Gibbs's definition. At low T, the curve is well fitted with 
ci{T /Tpc)~^ (the blue dashed line). The pion mass is extracted from the coefficient Ci of 
the 1/T fit at low temperatures. 



The pion mass in Eq. ()2-9p is shown in Fig. [10] as a function of T/Tpc- The results 
are shown for rrips/my = 0.6 and 0.8. In case of mps/my = 0.8, m-,^/T is well fitted 
with ci{T/Tpc)~^ at T/Tpc = 0.5 with ci = 4.060(6). The pion mass is extracted 
from the low temperature behavior as m.,^ = ciTpc, which is approximately 4Tpc. 
This large value is because of the present simulation setup. At high T, is almost 
linearly proportional to T, which is probably due to thermal effects. The decreasing 
"ips/w-v makes m-j^/T smaller. However, simulation data are not sufficient to obtain 

Cl. 

Hadron masses are extracted through the exponential damping in an usual 
method. In the present approach, the pion mass is extracted through 1/T behavior. 
It is important to note that the 1/T behavior is already obtained at Nt = 8. Hence, 
this approach may be an alternative method for hadron spectroscopy, which was 
discussed in detail in Ref.™'' 

§3. Low^ temperature limit of QCD 

In this section, we consider the behavior of the fermion determinant at low 
temperatures. For small ^, the fermion determinant becomes insensitive to /i as T 
decreases, and becomes almost independent of ^ at enough low T. We explain this 
behavior by using the properties of the eigenvalues of the reduced matrix. 

We first consider T- and /^-dependence of a reweighting factor. Then, we consider 
the low temperature limit with the use of the A'f-scaling law of the eigenvalues of 
the reduced matrix. We will derive two expressions for low temperature limit of 
the quark determinant; one is for small /i and the other is for large /i. The low 
density expression shows the /^-independence of the fermion determinant at T = 
for ^ < m,r/2. We will conclude that the /i-independence at T = for small is 
the consequence of the Nt scaling law of the eigenvalues of the reduced matrix. The 
other expression is its high-density counter part. We discuss the physical meaning 
and criterion for these limits. 
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We discuss the high density limit of the fermion determinant. In low temperature 
and high density limit, QCD approaches to a theory, where quarks interact in spatial 
directions with the ordinary Yang-Mills type of the gauge action. The corresponding 
partition function is invariant even in the presence of dynamical quarks. The 
fermion determinant becomes real and the theory is free from the sign problem. 

3.1. Fluctuation of the fermion determinant 

In this subsection, we consider the quark determinant at low temperature and 
finite fi. Figure [TT] is the scatter plot of lni?(^, /3o)(o,/3o) = iVj Indet Z\(/x)/ detZ\(0), 
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Fig. 11. The scatter plot of the reweighting factor on the complex plane. Left panel : hadron phase 
{T/Tpc = 0.93). Right panel : QGP phase(r/rpc = 1.08). The figures are taken from Ref.^^ 
The results are obtained for 8^ x 4 lattice with mps/mv ~ 0.8. 

which is so called reweighting factor. This shows how the quark determinant develops 
as iJ, is varied from a simulation point (/i = 0, /3o), where gauge configurations are 
generated. Here we fix a temperature (/3 = /3q). The left and right panels show the 
results for T/Tpc ~ 0.93(/3o = 1.8) and 1.08(1.9), respectively. The horizontal and 
vertical axes are the real and imaginary parts of Ini?, i.e., the exponent of \R\ and 
the phase of R. 
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Fig. 12. Left : the average of ln|/?| = A^/ In | det Zi(/x)/ det Zi(0) |, right : the average phase factor 
defined by (cos0), where 9 = arg(i?). The results are obtained for 8^ x 4 lattice with mps/mv = 
0.8. 
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The results are shown in Fig. [12] now as functions of T/Tpc for three values of 
fi. The magnitude of the reweighting factor is very small at low temperature, and 
rapidly increases near Tpc. As T decreases, the average of In \R\ approaches to zero at 
least up to na = 0.25. This behavior means | det ^(/i)/ det Z\(0)| ~ 1, and therefore 
the fermion determinant is insensitive to /i at low temperatures T/Tpc ~ 0.5. In the 
right panel, we plot the average phase factor {cos 9) with 9 = arg(i?). The average 
phase factor is close to one at high and low temperatures, and has minimum at 
slightly below Tpc. The sign problem is very severe in the vicinity of Tpc. At high 
temperatures, the average phase factor increases with increasing /x. On the other 
hand, the //-dependence of the average phase factor is small at low temperatures. 
The average phase factor approaches to one for the three values of fi. 
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Fig. 13. The quark chemical potential(/i) dependence of the average of In |7?| — Nf In | det A{jj,) / det A{0) \ . 



We found that the fermion determinant is insensitive to at low temperatures. 

However, the low-T behavior of the fermion determinant depends on //. Next, we 

consider the reweighting factor as a function of /i in Fig. [131 At high T(T > 0.76Tpc), 

the average value of In \R\ smoothly increases as fi goes to larger. A discontinuous 

change is found at fia = 0.5 ~ 0.6 at low T(T < 0.54Tpc). \R\ is approximately unity 

for small /i, and starts to increase at iia = 0.5 ~ 0.6, see the right panel. The onset 

of the /i-dependence of the fermion determinant is about /ia = 0.5 ~ 0.6(/i/r = 

4.0 ~ 4.8). Using the value of the pion mass obtained in the previous section, it 

corresponds to /x = 0.5m.,^ ~ O.Gm-^. This is consistent with a well known result 

that the fermion determinant is independent of for /x < at T = 0, see e.g. 

pj^gfg 23),26),27.),68J 

Phenomenologically, it is expected that the /x-independence at T = continues 
up to /X = Mjv/3, where M^r is the mass of the nucleon. This problem was raised 
long time ago in the Glasgow method, see e.g. Ref.— ' There would be several 
reasons for the discrepancy between the critical value of /x and Mjsf/3. For instance, 
the importance sampling at /x = may cause this discrepancy, because the quark 
chemical potential is equivalent to the isospin chemical potential at /x = 0. 

The deviation of the reweighting factor, cj^, is shown in Fig. [TH The deviation 
of the reweighting factor reaches the maximum near the crossover transition point 
Tpc both in the magnitude (left panel) and the phase (right panel), and decreases 
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Fig. 14. The fluctuation of the reweighting factor as a function of T/Tpc- = ^ ~ a;)^,a; = 

^^Hi^M where a:: = Re[ln /3o)(o,0o)l f'^'^ ^^^^ panel and x = Im[ln /3o)(o,^o)] ^'^'^ right 
panel. The results for Nt — are taken from 



as the temperature is away from Tpc. The peak becomes prominent as ^ increases. 
In Fig. nn we showed the reweighting factor for low T{T/T^c = 0.93) and high 
r(1.08). In the vicinity of Tpc, gauge configurations visit low-T states and high-T 
states, which results in the peak of the fluctuation. The real part of In \R\ reaches the 
maximum at Tpc, while the imaginary part reaches the maximum at slightly below 
Tpc. Hence the sign problem is most severe at slightly below Tpc, see also the right 
panel of Fig. [T2j 



Fluctuation of IVIagnitude Fluctuation of Phase 




Fig. 15. The fluctuation of the quark determinant as a function of T/Tpc and the Gaussian fit. See 
also Fig. [H 

In Fig. [T5l we focus on the low temperature region. As T decreases, the de- 
viation rapidly decreases for both real and imaginary parts. We find that the low 
temperature behavior is well fitted with the Gaussian function, which is shown in 
the solid lines in Fig. [T5j This implies that fluctuation of the quark determinant is 
exponentially suppressed as the temperature decreases for small //. 

We have seen that the fermion determinant is insensitive to /i at low temper- 
atures for jj, < m7r/2. As we have discussed in § 12.51 the eigenvalues follow the 
A^t-scaling law. Let us consider an large eigenvalue A (|A| > 1) and describe its coun- 
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terpart as 1/A*. As increasing Nt or decreasing T, the large eigenvalue A becomes 
lager and the smaller eigenvalue 1/A* smaller. If /i is fixed at a small value, is 
0(1). This leads to the scale separation |1/A*| ^ ^ ^ |A|. The contribution of each 
eigenvalue to the fermion determinant is approximated as ^ + A ~ A and ^ + 1/A* ~ ^. 
This causes the //-independence of the fermion determinant at low temperatures. We 
will turn back to this point in the next subsection. 

The average phase factor {cosO) is close to one at low temperatures as well as 
high temperatures. Does it mean the feasibility of MC simulations at low T ? At high 
T, although the fluctuation is small, the fi dependence of the fermion determinant 
is non vanishing, while the //-dependence almost disappears at low T. It is already 
known that the neglecting phase leads to the phase transition at /i = mT^/2. The 
phase of the determinant plays an important role to go beyond ■m-^j'l. The careful 
analysis would be required to evaluate the phase of the determinant, which may 
cause another difficulty at low temperature simulations. 
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Fig. 16. The fluctuation of the quark determinant on 10 x 4. See also the results for A^'s = 8 in 
Fig. 1141 G is given in the caption of Fig. 1141 



Note that the above results are obtained for fixed spatial lattice size Ng. As we 
have mentioned, the fermion determinant depends on A^^, although the eigenvalues 
of the reduced matrix are insensitive to Ng. Finally in this subsection, we consider 
the volume dependence of the reweighting factor. The result for Ng = 10 is shown 
in Fig. [T6j See also the result for A^, = 8 in Fig. [TH Increasing A^, from 8 to 
10, the maximum value of the deviation is almost twice both in the left and right 
panels, which is approximately equal to the volume ratio Vs{Ns = lG)/Vs{Ns = 8) = 
lOVS^ ~ 2. 

A question arises on the phase of the fermion determinant at the low temper- 
ature limit Nt ^ oo and thermodynamical limit A^g — t- oo. As we have shown, 
decreasing temperature suppresses the phase, while increasing the spatial volume 
enhances the phase. What happens for the phase at the low temperature limit and 
thermodynamical limit ? The result depends in general on the order of taking these 
two limits. In this paper, we gave only a result in Fig. [12] which was obtained for 
a fixed lattice volume. We will extend the present study for various sets of lattice 
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sizes Ng and Nt in the next step 0- 

3.2. Low temperature limit of the fermion determinant 

According to the A^^-scaUng law, we parameterize the larger half of the eigenval- 
ues as A„ = Z^'e*^". Using the pair nature, the smaller half of the eigenvalues are 
described as 1/A* = /~^'e*^". Now, the quark determinant is parameterize as 

Afrcd/2 Af,cd/2 

detZ^ = Cor^-^/' n (^ + ^n'e''") n (e + C^'e^^")' (3-1) 

n=l n=l 

where we set c± = 1 for simplicity. In Eq. (j3-ip . we divide the product into two 
parts and replace the smaller eigenvalues with 1/A*. 

In the low temperature limit T = l/{aNt) — > 0, (A^ — > oo), large and small 
eigenvalues follow 

e^ocGO-'^^^O. (3-2) 

Since the fugacity ^ = exp(— /ioA^) also decreases with the same exponent as the 
small eigenvalues, det in the low-T limit depends on /i. 

I) Fixed /i/T. In this case, the fugacity ^ is constant, and the smaller eigenvalues 
decrease faster than ^ ^ Then, the quark determinant is reduced to 

Afred/2 

del A = Co n (3-3) 

n=l 

where the product is taken over the larger half of the eigenvalues [ A„ | > 1. det A[n) 
is independent of /i, and therefore there is no sign problem in this limit. 

II) Fixed /ua. In this case, the fugacity decreases with the same exponent as the 
smaller eigenvalues. This case is further classified into two cases corresponding to 
the magnitude relation of exp(^a) and In- We introduce a typical magnitude /, 
which is used to compare the eigenvalues and exp(/ia). This may be the average 
of /„ or smallest one of depending on the distribution of eigenvalues. We will 
turn back to this point later. 

a) exp(/ia) < /. In this case, the smaller eigenvalues decrease faster than the 
fugacity; ^ » (/n)~^*- We obtain 

A^rcd/2 

detZ\ = Co n (3-4) 

n=l 

This is equal to Eq. (j3-3p . which implies that the quark determinant remains 
unchanged up to /io < In/ in the low temperature limit. Namely, detZ\(/i) is 
independent of /x at low temperature and small ^ regions. This is nothing but 
the Silver Blaze phenomena discussed above. 

*•* SplittorfT and Verbaarschot investigated the average phase factor at low temperature in chiral 
perturbation theory,!^ and discussed this problem. 
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This is analogous to a situation in Fermi statistics. If T is smaller than a lowest 
excitation energy of a system, then the inclusion of small fi can excite no quark 
in excited energy levels, the system remains to stay at the lowest energy state. 
This can lead the system to be independent of /i. We discuss this point in Ref.!^ 

b) exp(/ia) > /. In this case, the fugacity decreases faster than the smaller eigen- 
values; ^ < Then, we obtain 

detZ\ = r^-d/2(^^^g^Q 

Nt 

= r^-/2[]det(A) 

i=l 

Nt 

= -Q^g^(^ac,M^(^^^^^.) ^au _ r^^Six-y)). (3-5) 

■t=l 

where we first use Eq. (|2-5|) . then substitute Eq. (|2-4c|) . In the last line, we use 
detC/ = 1. In Eq. ()3-5p . B{ti) contains spatial hops in the i-th time slice, but 
does not contain any temporal hopping terms. The jj, dependence comes only 
from the overall factor ^~^rGd/2 _ Qx.p{2NcN^ ^/T). This is the highest order 
term in the fugacity expansion and means that all the states are occupied by 
quarks. 

As we have discussed, det Q = 1 and Cq is real, therefore Eq. ()3-5p is real and free 
from the sign problem. The fermion determinant is real in the low temperature 
limit both for large and small fi. However, the fermion determinant is given 
by the different expressions at large and small n, which may suggest that two 
different states exist at T = 0. 

One may think to determine the critical value of ^oa = In/. However the deter- 
mination of / is nontrivial task because of the finite width of the distribution of the 
eigenvalues. If we employ the eigenvalue closest to 1 for I, then T"^* = m.a'K\x,^\<:i |A„|. 
Using Eq. ()2-9p . we obtain 

In I = ij^a = am-jr/ 2. (3-6) 

The fermion determinant is independent of /i for < m^/2. This value of / corre- 
sponds to the phase transition point for a pion condensation observed in the phase 
quench simulations. If we employ the average value of the larger half of A„ for /, we 
approximately obtain I ^ (200)^/^. We obtain fia ~ Inl = ;|ln200 ~ 1.32 which is 
much larger than and beyond the lattice cutoff //a = 1. 

The criterion / may be different for the low density limit and high density limit, 
depending on the eigenvalue distribution at T = 0. The situation of the low tem- 
perature limit is shown in a schematic figure Fig. [T71 If fi is smaller than m.j^/2, the 
fugacity is located within the gap. Taking Nt ^ oo leads to |1/A*[ ^ ^ ^ |A|, (for 
|A| > 1), which corresponds to (Il-a). The fermion determinant is independent of fi 
in this case. Increasing ^, the fugacity becomes comparable to the smaller half of 
the eigenvalues, which causes the fj, dependence. If /i goes beyond a certain value. 
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high density 
limit 



-/ifl=(ln^)/]V, 




A\/N, 



Fig. 17. Schematic figure for the low temperature hmit. The sohd fine denotes the histogram 
of \n\\\/Nt, see also Fig. [H] Dotted vertical lines denote the behavior of —fi/T = In^ with 
increasing /x. 



which is probably given by the minimum eigenvalue, then the low temperature limit 
(Il-b) is obtained. The fermion determinant has a trivial ;U-dependence in this case. 
According to the above discussion, the criterion would be given by 

> max|A|, for (Il-a), (3-7) 

A|<1 

C<min[A[, for (Il-b). (3-8) 

The criterion for (ITa) is related to the pion mass. The criterion for (Il-b) may also 
have a similar interpretation. As we have discussed in § 12.21 the eigenvalues of the 
reduced matrix is related to the free energy of quarks. According to the discussion 
there, the minimum eigenvalue is related to the highest energy state of a quark. 
Hence, the low temperature limit (ITb) is obtained if fi is sufficiently larger than the 
highest energy state of a quark. 

A question arises if the high density limit reflects a real physics of QCD or just a 
consequence of lattice artifacts, because the highest energy state probably depends 
on the lattice spacing a and is a consequence of UV-cutoff. The a-dependence of 
the minimum eigenvalue is understood from the left panel of Fig. SI The increase 
of T means the decrease of a, which follows from T = {aNt)~^. It turns out that 
the maximum eigenvalue Amax becomes larger as a decreases. This means that the 
minimum eigenvalue Amin becomes smaller with decreasing a because of the relation 
Amin = l/'^max- The minimum eigenvalue is described as Amin ~ e~^/-^ assuming the 
eigenvalues correspond to free energies of a quark. With decreasing a, Amin becomes 
smaller, which means F becomes larger. Thus, the highest energy state depends on 
the lattice spacing a. Investigations of the eigenvalue distribution for finer lattices 
would clarify if the limit (H-b) corresponds to the low temperature and high density 
limit of QCD. 

It is also important to consider the dependence of Amax and Amin on the quark 
mass and lattice volume as well as the lattice spacing. In the present work, the value 
of Amax depends on the quark mass (Fig. [6]), while it is insensitive to the lattice 
volume A'^^ (Fig. E]). 
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3.3. QCD at low temperature and high density 

In this subsection, we study further the low temperature and high density Umit 
(Il-b). Even though the hmit (Il-b) may be the consequence of the lattice cutoff, it 
is meaningful to consider this case. For instance, it can be used to generate gauge 
configurations at high density regions. So far, there is no case where direct MC 
simulations are feasible and the quark number density is high. If gauge configurations 
are generated at high density regions by direct MC simulations, they may provide 
valuable information, e.g., for multi-ensemble reweighting. 

Using Eq. (j3-5p . we obtain the low-temperature and high-density limit of the 
QCD partition function 

lim Zacil^, T) = e^^f^^^st^/^ f VU (det Z\(^)|t-,o)^^ 6"^^, (3-9a) 

T^O J 
Nt 

detA{fi,T)\T^o = l[det {B^'^^%x,y,t,) - 2k r''_'6{x-y)) , (3-9b) 
j=i 

where the definitions of B, r± etc were given in § 12.11 

In Eq. (|3-9p . the gluon part remains unchanged and is given by the ordinary 
Yang- Mills action, while the fermionic part is different from the ordinary QCD action. 

Quarks interact only in spatial directions, where no interaction exists in the 
temporal direction. 

Equation (|3-9bp is also different from the naive spatial fermion matrix B, but 
contains a constant term 2k with the projection operator r± = (r ±74)72. The quark 
chemical potential appears only in the bulk factor exp{2NfNcNgfi/T), which gives 
the quark number density, 

(n) = 2NfNc, (lattice unit) (3- 10a) 

(X 2NfNcH^, (physical unit). (3-lOb) 

Since Eq. (|3-9bp contains no temporal hopping term, the partition function is 
symmetric even in the presence of dynamical quarks. Naively, it is expected that 
a deconfinement transition occurs if baryon number density is large enough to cause 
the overlap of baryon's wave functions, where effective degrees of freedom would be 
quarks rather than hadrons. If this is the case, symmetry would be an exact 
symmetry of QCD in extremely high-dense matter. 

Another high density limit was proposed in Ref.'^S with an approximation that 
the quark mass and chemical potential are simultaneously made large. In the ap- 
proximation, the quark mass depends on the chemical potential. In the present case, 
the fine tune of the quark mass is not needed in taking the low-T and high-/i limit. 

Equation ()3-9p is realized at extremely large fi. It would be very difficult to 
access such a high density region in experiments. Nevertheless there are several 
interests in the high density limit (Il-b). Theoretically, it is interesting to consider 
the nature of QCD at high density regions : confinement, chiral symmetry and 
color superconductivity. In lattice QCD simulations, the knowledge on important 
configurations at high density regions would be valuable information. For instance. 
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such configurations can be used in multi-ensemble reweighting method, or they may 
be used for a reweighting from high density regions in order to find a QCD phase 
transition line at low temperatures. 

It is important to consider the numerical feasibility. As we have mentioned, 
whole the determinant is real, and therefore sign free. Although it is free from the 
sign problem, it needs large Nt to take the low temperature limit, which requires 
a large numerical cost. This increase of the computational time may be suppressed 
by using the property of the fermion determinant in the low temperature limit. In 
the low-T and high-^ limit, the fermion determinant is expressed as the product of 
the Nf block determinants. If each block determinant is real, it may be possible to 
evaluate the block determinant by the Gaussian integral with the pseudo-fermion 
field with smaller dimension. Instead, the Gaussian integral appears Nt times. The 
real positivity of the block determinant is necessary to implement this idea. We leave 
the proof of this expectation in future works. 

§4. The partition function on the complex plane 

The determination of the confinement /deconfinement phase boundary is an im- 
portant issue in the study of the QCD phase diagram. A canonical formalism and 
Lee- Yang zero analysis are useful approaches to identify a phase transition point, and 
have been investigated in Refs.[33.ED,llD,[I3,ll3,C3,C3,E3,[Z3,[Z3,[Z3,B3,ll3,Cl How- 

ever, some difficulties were pointed out.!^'!^ For instance, the sign problem causes 
a fictitious signal in Lee- Yang zero analysis, which makes it difficult to distinguish a 
physical phase transition and fictitious signal.'^ 

In this section, we consider the canonical formalism and Lee- Yang zero theorem, 
where the reduction formula is useful. The purpose of this section is to propose a 
method to identify a /x-induced phase transition by using a temperature dependence 
of canonical partition functions and of Lee- Yang zero trajectories. The method 
provides a qualitative way to distinguish a crossover and first order phase transition, 
and would be useful in case that the sign problem causes a difficulty in ordinary 
methods such as finite size scaling of the Lee- Yang zero near a positive real axis. 

Canonical partition functions Zn with the quark number n are obtained from the 
fugacity expansion form of the fermion determinant. Then, we consider the temper- 
ature dependence of two quantities: the canonical distribution Z„ and trajectory of 
Lee- Yang zeros. They show characteristic changes as the temperature decreases from 
high T to low T. We discuss the relation between their temperature dependences 
and the existence of a /^-induced phase transition. 

In § 14.11 we give a brief overview of the canonical formalism and Lee- Yang zero 
analysis. Fugacity coefficients of the fermion determinant are calculated in § 14.21 
Canonical partition functions are calculated in § 14. 3t where we employ the Glasgow 
method.l^ Lee- Yang zeros are calculated in § 14.41 
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4.1. Brief overview 

According to statistical mechanics, a grand canonical partition function is de- 
scribed as a superposition of canonical partition functions, 

N 

Z{T,f,) = ^Zn{T)C, (4-1) 

n=l 

where Zn describes a canonical partition function with a fixed particle number n, 
and N is the maximum number of the particle. In thermodynamical limit, N — )• oo. 

Several methods are available for the determination of phase transition points. 
For instance, a phase transition point is identified by the convergence radius of 
Eq. (|4-ip . or the finite size scaling analysis of a Lee- Yang zero Z{iJ,) = 0, the Maxwell 
construction for an S-shape in a ii-n diagram, etc. 

The grand partition function Eq. (|4Tp can be considered as a polynomial of ^ 
at a given temperature. Phase transitions result from discontinuities of derivatives 
of the free energy. In general, Z{^) contains the A^-roots on the complex ^ plane. 
Since the canonical partition function Zn is real positive for all n, no root exists 
for real positive values of ^. Hence, all the roots are distributed somewhere on the 
complex ^ plane except for the positive real axis. Lee and Yang showed ESJ-EO) that 
in the case a phase transition occurs, zeros of the grand partition function approach 
to the positive real axis in the thermodynamical limit V oo. Thus, the phase 
transition is described by zeros of the grand partition function, which are called the 
Lee- Yang zeros. The order of the phase transition is distinguished by considering the 
trajectory formed by zeros near the positive real axis.l^ In Ref.'^ the application to 
non-equilibrium systems was also discussed. Stephanov investigated the properties 
of Lee- Yang zero of QCD by using the scaling and universality.'^ In lattice QCD 
simulations, the thermodynamical singularities of Lee- Yang zeros were investigated 
in e.g. Refs.[21'ESj),^,^,M) 

4.2. Fugacity expansion of the fermion determinant 
We start from the fugacity expansion of det A{p), 

Afrcd/2 

det A{fi) = Co Yl (4-2) 

n=-Afred/2 

which is defined in Eq. ()2-6bp . 

Before seeing numerical results, we discuss the numerical procedure for the calcu- 
lation of c„. There are two difficulties in the determination of c„; numerical precision 
and overflow/underflow due to the applicable range of the double precision. 

First, the expansion of Eq. ()2-6ap requires an enormous amount of calculations 
because of the larger number Ared, which easily causes a loss of precision. A recursive 
algorithm with a recurrence relation was used to expand Eq. (|2-6ap . A Vandermonde 
matrix approach was not effective because of the large dimension and of the existence 
of close eigenvalues. Divide and conquer algorithms are best in respect to the num- 
ber of the calculation steps. Second, Cn span a wide range of magnitude, from order 
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0(1) to 0(e^^'=^'' ). A huge number such as 0(e^^'=^» ) exceeds the maximum value of 
double precision, where an ordinary double precision variable is not applicable. Nu- 
merical libraries such as FMLitP^ are available for this calculation. Other libraries 
have been also used in literature, see e.g., Ref.l^ However, the use of libraries may 
become an obstructive factor of fast computation. Therefore, we developed a new 
variable.SS' These procedures to calculate c„ are explained in Appendix. O 




n 



Fig. 18. The magnitude of all the fugacity coefficients for various temperature. The lattice size is 
8'^ X 4. Data is taken from one configuration for each temperature. 
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n n n 

Fig. 19. The argument of the fugacity coefficients for various temperature. Data are shown for 
small quark number sector. The lattice size is 8^ x 4. Data is taken from one configuration. 

Now we proceed to the numerical results of c„. Figures [18] and [19] show the 
modulus and argument of the fugacity coefficients c„ for a gauge configuration. The 
simulation setup was given in the section [2.31 The magnitude |c„| spans over the 
wide range of order from 0(1) to 0{e^^ ). At /i = 0, detZ\(/x) is dominated by 
several c„ near n = 0. The argument of Cn shows complicated n dependence, as 
shown in Fig. [TD] where arg(c„) is defined for —it < arg(c„) < vr. Qualitatively, 
arg(cn) tends to oscillate with higher frequency as the temperature becomes lower. 
Calculating c„ for 400 configurations, we found that \cn\ was stable for the change 
of configurations. On the other hand, arg(c„) rapidly changes for configuration by 
configuration, which leads to the cancellation of Cn in the ensemble average. This 
cancellation becomes more severe at low temperatures. We will see this point in the 
next subsection. 

The fugacity coefficients c„, satisfy the relation cl„ = c„, as a consequence of 
the 75 hermiticity. Hence, positive and negative n describe the winding number of a 
quark and an anti-quark around the temporal cylinder, respectively. This is realized 
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as the reflection symmetry |c_„| = |c„| with respect to n = for the absolute value, 
which is well satisfied see Fig. [THJ The relation for the phase, which is given by 
arg(c_„,) = — arg(c„), was numerically satisfied for several hundreds c„ near n ~ 
and near n ibA'red/2. For instance, in the left panel of Fig. [191 arg(c_n) = — arg(c„) 
for n = ~ 300, while arg(c_n) / — arg(c„) for n ~ 500. Fortunately, the quark 
determinant is dominated by coefficients near n = for = 0. 

4.3. Canonical partition function 

The grand canonical partition function of QCD with Nf flavors can be also 
expanded in powers of the fugacity, 

ZGc{f^,T)= Zc{n,T)C, (4-3) 



where Ng = NfNj-ed/^ = 2NfNcN^ is the maximum quark number which can be put 
on the lattice. Zc(n,T) is a canonical partition function with a fixed quark num- 
ber n. Using the Fourier transformation,S3'nSli'H3 the canonical partition function 
is obtained 

Zc{n) = r^"" d<Pe-'^^ZGc{l^ = i/x,), (4-4) 

where cp = l^-i/T- We have used the Roberge- Weiss periodicity. 

In this work, we employ the Glasgow method, which is based on the reweighting 
in fi and reduction formula, for the calculation of the canonical partition functions. 
It was pointed out that it suffers from the overlap problem.^ In this work, we focus 
on the properties of the canonical partition functions and Lee- Yang zeros as a first 
step, and leave the improvement of the overlap for future works. 

Substituting the reduction formula into Eq. (|4-4|) . the canonical partition func- 
tion is given by 

Here dn are the fugacity coefficients of the two-flavor determinant, and (•)o denotes 
an ensemble average for gauge conflgurations generated at fi = 0. dn is determined 

by using the recursive algorithm for (r^^^''^^ 11^=^ + A„)) ^ = Enio'"" ^nC ■ 
Although it is possible to obtain dn in terms of c„, this method was slower than the 
above procedure. 

The canonical partition functions Zn are shown up to \n\ =30 for three tem- 
peratures in Fig. [20l Z„ must be real positive, i.e., Im[Z„] = 0. The results for 
high temperatures satisfy this condition. Even at low temperature, Im[Z„] is zero 
within errorbars for most n. However, the errorbars are large, which is caused by 
the fluctuation of the phase of the fugacity coefficients c„ for configuration by con- 
figuration, which is the sign problem in canonical approaches.!^ The signal-to-noise 
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Fig. 20. The canonical partition functions Z„ for three temperatures. Left : real part, right: imag- 
inary part. In the left panel, the solid lines are exp(— ai|n|) for T/Tpc = 0.93 and exp(--a_f/n^) 
for T/Tpc = 1 and 1.35. The results are obtained by using the Glasgow method with the use of 
400 sets of the eigenvalues of the reduced matrix. 

ratio for Re[Z„] becomes small for large n. This is probably because of the impor- 
tance sampling at ^ = 0. The average value of the quark number density is zero at 
fi = 0, and therefore configurations generated at /i = have less overlap with large 
quark number sectors. The canonical partition functions are obtained up to about 
|?7-| = 30 in the present simulation setup. The calculation of Z„ for larger n needs an 
improvement of overlap or increase of statistics. However, it should be noted that 
the real part of the canonical partition functions show a characteristic temperature 
dependence even for |n| < 30. 

The canonical partition function Re[Z„] exponentially decreases as |n| increases 
for all the three temperatures. The width of the distribution of Re[Z„] becomes broad 
at high temperature, which is consistent with the increase of the effective degrees 
of freedom at high T. The n-dependence of Re[Z„] changes at high temperatures 
{T/Tpc = 1,1.35) and low temperature {T/Tpc = 0.93). Qualitatively, Re[Z„] is 
approximated by a Gaussian function exp{—aHn'^) for T > Tpc, and to a function 
exp(— Oilnl) for T < Tpc, where aH,L are parameters. To be precise, the result 
agrees with the Gaussian function for T/Tpc = 1-35 up to |n| = 30, while there is a 
deviation from the Gaussian function for jnj = 20 ~ 30 for T/Tpc = 1- There is also 
deviation from the function exp(— aL|n|) for T/Tpc = 0.93. In case of exp(— a^lnl), 
the partition function becomes a geometric series. 

Triality nonzero terms mod(n, 3) 7^ do not vanish, which is a consequence of 
the importance sampling at /i = 0. Those terms can be eliminated by using the 
Roberge- Weiss periodicity.— ' We discuss this point in Appendix. |Bj The triality 
nonzero terms change neither the trajectory of Lee- Yang zeros nor the convergence 
radius of the fugacity polynomial as long as Z„ with zero- and nonzero-triality is 
described by the same function of n. Triality nonzero terms change the density of 
Lee- Yang zeros on its trajectory. However, this effect vanishes if thermodynamical 
limit is taken. Hence, we keep the triality nonzero terms, which does not cause 
any problem in the following discussion. Here it is interesting to note that the low- 
temperature high density limit of QCD may suggest the importance of the other Z3 
sectors. If all the Z3 sectors are visited in MC simulations, the triality nonzero terms 
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would disappear. 

The n-dependence of the canonical partition functions change from high tem- 
peratures to low temperatures. Now, we discuss the relation between this change 
of the n-dependence of Re[Z„] and a finite density phase transition. We found that 
Re[Z„] is approximately given by 

( p-ah\n\ (T < T ) 

Re\Zc{n)] - T . ' (4-6) 

Assuming that these distributions hold for larger n and that the imaginary part is 
sufficiently small, the convergence radius r of Eq. (j4-3p is given by 



r ^ = lim 



Z{n + 1) 



Z{n) 



(T > Tpc). 



(4-7) 



The Gaussian function for T > Tpc suggests no phase transition, while the function 
Q-a.\n\ £qj. rp ^ j,^^ suggests the existence of a phase transition at finite fi. Note 
that the result for T = Tpc shows the Gaussian behavior, which is consistent with 
the absence of a /i-induced phase transition at T = Tpc. Thus, the shape change of 
the canonical distribution is related to exisitence or absence of a /x-induced phase 
transition. 

The deviation from Eq. ()4-6p observed for T = Tpc and T = 0.93Tpc may be 
important in the study of CEP. The shape of the canonical distribution including 
the large-n part and the deviation from Eq. ()4-6p can contribute to higher-order 
moments, such as skewness and kurtosis. In the present analysis, we employed the 
one-parameter reweighting, where the overlap problem becomes severe for large n 
parts of Zn- The improvement of the overlap is needed in order to apply the above 
discussion to cases where the tail part of the canonical distribution is important. 

4.4. Lee-Yang zeros 

In this subsection, we consider the Lee- Yang zeros. Several methods are available 
for the calculation of Lee- Yang zeros. For instance, it is possible to search zeros of the 
left hand side of Eq. (|4-3p . It is also possible to solve roots of the fugacity polynomial, 
which is the right hand side of Eq. ()4-3p . Here, we adopt the latter approach by using 
the canonical partition functions obtained in the previous subsection. 

In order to calculate roots of the fugacity polynomial. We consider the trunca- 
tion of the fugacity polynomial in order to calculate roots of the right hand side of 
Eq. (|4-3p . because the order of the fugacity polynomial Ng is large. In general, it 
is not allowed to truncate the polynomial in the vicinity of phase transition points. 
However, the truncated polynomial can reproduce the trajectory of Lee- Yang zeros 
in the case of the geometric series. First, we discuss this point. 

We divide the fugacity polynomial into small n and large n parts 

zgM= zc{n)e+ (4-8) 

\n\<M \n\>M 

where M is an integer. Although Re[Z„] rapidly decreases as |n| increases, the sum 
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of the higher-order terms significantly contribute to Zgc{i^) in the vicinity of the 
convergence radius, i.e., phase transition points. 

At high r, the higher-order terms do not affect the convergence property because 
of the Gaussian shape. At low T, the e~°l"'l-shape provides the nonzero convergence 
radius, where the sum of higher-order terms is significant near transition points. 
However, if Re[Z„] = exp(— a|n|) even for large n, the trajectory formed by Lee- Yang 
zeros does not depend on the maximum order of the fugacity polynomial. This is a 
consequence of the geometric series. For instance, considering l+x+x'^+- ■ ■+x^ = 0, 
its roots lie on the unit circle. The order changes the density of the roots on the 
circle, but does not change the trajectory formed by roots. 

In fact, the fugacity polynomial with Re[Z„] = exp(— a|n|) is different from a 
naive geometric series, because it contains both the quark and anti-quark compo- 
nents. However, as we will show below, the trajectory of the roots for the quark and 
anti-quark sectors are separated from inside and outside of |^| = 1 because of the 
symmetry between quarks and anti-quarks. 



E 




-3-2-10123 

Fig. 21. The distribution of Lee- Yang zeros, which are roots of the fugacity Polynomial shown in 
Fig |20l The roots are obtained for the leading order terms for M = 24. Solid lines are the unit 
circle. 

Accordingly, the Lee- Yang zero trajectory can be reproduced via the truncated 
polynomial on the assumption on Z„. Now, we consider the numerical result. We 
employ the data of Z„ up to M = 24, because the signal-to-noise ratio becomes large 
for n > M at T/T^c = 0.93. Although this value of M is small, it corresponds to 
the baryon density ~ 2 [fm^^]. IMSL Library was used for the calculation of the 
roots of the fugacity polynomial. The result is shown in Fig. [2TJ Corresponding 
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to the change of Zn, the Lee- Yang zero trajectory also changes its shape from high 
temperatures to low temperatures. The zeros are approximately distributed on two 
circles at T/T^c = 0.93 and one circle with two branches at T/T^^ = 1 and 1.35. 
The trajectory is symmetry with regard to the unit circle because of the charge 
conjugation symmetry between the quark and anti-quark. 

The results for the Lee- Yang zero trajectories were obtained by using the data 
of Zn shown in Fig. [201 where both the real and imaginary parts were considered. 
We did not use Eq. (j4-6p to obtain Fig. [2TJ Hence, the trajectories in Fig. 1211 contain 
the deviation of Z^ from Eq. ()4-6p . The two-circle trajectory at low T is similar to 
a typical behavior of geometric series with positive and negative n components. If 
the polynomial is an ordinary geometric series with positive n components, then the 
roots lie on a circle. Adding the negative n components with the charge conjugation 
symmetry, then the roots locate on two circle. Similarly, the two-branch trajectory is 
a typical behavior of the Gaussian distribution. Therefore, the trajectories obtained 
from the data of Zn qualitatively agree with the trajectories obtained from Eq. ()4-6p . 
which suggests that Eq. (j4-6p is a good approximation for Z„ at least for small n. 
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Fig. 22. The distribution of Lee- Yang zeros on the complex /i/T plane. 



In Fig. [22| we have shown the Lee- Yang zero distribution for T/Tpc = 0.84 and 
0.93 in the complex /i/T plane. The results are invariant under fi/T —^/T. 

In the Lee- Yang zero theorem, the phase transition point is obtained from the 
finite size scaling analysis of the zero nearest the positive real axis. In the present 
approach, we have made the assumption on Z„. If the assumption is valid, then 
Lee- Yang zeros are on the same trajectories. Then, the zeros would approach to the 
positive real axis as the volume increases in case of a phase transition. The phase 
transition point can be estimated in the canonical formalism, by using the Maxwell 
construction for the S-shape of the /U-n diagram .'^3' 'S3 

Following the assumption Eq. ()4-6p . we estimate a phase transition point. Here 
we should note that the overlap problem causes a loss of reliability. From the linear 
fit Re [/i/T] = const., we estimate the location of the intercept of the trajectory and 
real positive axis, and obtain 



^c/T 



0.97(3) 
0.70(2) 



r/Tpe 
T/T^c 



0.84 
0.93 



(4-9) 
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Fig. 23. A conjectured phase boundary for the first order phase transition. The data are obtained 
from Fig. (|22p as the intersection of the linear fit and positive real axis with the assumption on 



The result is mapped onto the QCD phase diagram in Fig. [23j Here the errorbars for 
Hc/T were obtained from fit- We also estimated the errorbars for T/Tpc, which is 
taken from.ES The result is almost consistent with previous studies with staggered 
fermion^iS at fi/T = 0.7, and undershoots those results at /i/T = 0.97. 

Here, it is important to consider the applicable limit of the Glasgow method. 
This can be done by using the imaginary chemical potential. The pure imaginary 
chemical potential region is located on the unit circle on the complex fugacity plane. 
The Roberge- Weiss (RW) endpoint should be located on Im[/i/T] = -k/3. It turns 
out that the singularity on the unit circle appears about Im[/x/T] ~ 7i"/4, which 
is smaller than the value expected from the RW endpoint. This would imply that 
the overlap problem becomes severe for |^/T[ > 7r/4. The phase transition points 
in Fig. [23] are almost corresponds to this limit. In order to determine the phase 
transition point, the configurations should be improved to obtain a better overlap. 

In this section, we have presented an approach for the study of the QCD phase 
boundary. Although our analysis is at fundamental one, we found that the canoni- 
cal distribution and Lee- Yang zero trajectory distinguish the crossover behavior at 
Tpc and first order behavior at low T. It is useful to consider the canonical parti- 
tion function together with standard techniques for the phase transition, which may 
provides complementary information to identify the phase transition point. 

The assumption used in the present analysis should be examined further. In 
particular, the improvement of the overlap and the finite size effect on Z„ are im- 
portant task. We showed the results for Ng = 10 in Fig. [211 The increasing Ng 
causes the broadening of the canonical partition functions. Although the eigenvalue 
distribution of the reduced matrix is insensitive to Ng, canonical partition functions 
are sensitive to Ng. Since Nq for Ng = 10 is twice as larger as that for Ng = 8, 
we employed M = 48 for Ng = 10. The trajectory of the Lee- Yang zero was not 
affected largely by the increase of Ng. However, the lattice volume is still small for 
Ng = 10, and the finite size effect may appear for larger lattices. The finite size effect 
should be investigated for larger lattices. Note that in Fig. [231 the result for Ng = 8 
and 10 were obtained with the same number of statistics. Since the sign problem 
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Fig. 24. Volume dependence of the canonical partition functions (left panel) and Lee- Yang zero 
distribution (right panel). 

becomes more severe for lager lattices, it is also important to increase the number 
of the statistics to study the finite size effect. 

§5. Summary 

We have studied QCD at nonzero chemical potential and temperature in the 
lattice QCD simulations. We particularly focused on the low temperature regions of 
the QCD phase diagram, and studied several issues with the use of the dimensional 
reduction formula of the fermion determinant. 

In § [21 we studied the reduction formula of the Wilson fermion determinant and 
showed several properties of the reduced matrix. The reduced matrix is interpreted 
as the transfer matrices or the generalized Polyakov line, and the eigenvalues of 
the reduced matrix is related to the free energy of dynamical quarks. The angular 
distribution of the eigenvalues manifests the or confinement properties of QCD. 
The eigenvalues form the gap, which is related to the pion mass and therefore chiral 
symmetry breaking. We found an indication that the eigenvalues of the reduced 
matrix is scaled by the temporal size as A ~ for |A| > 1 and A ~ l^^^ for |A| < 1. 
The Nf scaling law controls the temperature dependence of the fermion determinant, 
and therefore is the important finding. 

In section [3l we studied the property of the fermion determinant at low T and 
finite fi. We showed from the lattice simulations that at T/Tpc ~ 0.5 and for /j, < 
771.^-/2, the fermion determinant is insensitive to and the average phase factor 
(cos 9) approaches to one. The fluctuation of the reweighting factor, the ratio of 
the determinant at /i = and /i 7^ 0, exponentially decreases as the temperature 
decreases. Using the A'^-scaling law, we also derived the //-independence of the 
fermion determinant at T = for fi < m7r/2. Hence we concluded that the fermion 
determinant is /^-independence at low T and for /i < m7r/2, which is the consequence 
of the A't-scaling law of the eigenvalues of the reduced matrix. 

Extending the low temperature studies further, we have considered the low- 
temperature limit of the fermion determinant and of QCD, and obtained two ex- 
pressions of the quark determinant; one is for low density and the other is for high 
density. Low density expression corresponds to //-independence for fj. < m-,^/2. The 
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other expression is its high-density counter part. 

We discussed the nature of the high density and low temperature Hmit of the 
QCD partition function. In the case, QCD approaches to a theory, where quarks 
interacts only in the spatial direction with the ordinary Yang-Mills type of the gauge 
action. The corresponding partition function is Z3 invariant even in the presence 
of dynamical quarks. Furthermore, the fermion determinant becomes real and the 
theory is free from the sign problem. 

In § m we studied the canonical formalism and Lee- Yang zero theorem. We 
have shown that the canonical distribution and Lee- Yang zero trajectory show char- 
acteristic changes when T is varied from high to low temperatures. The canonical 
distribution is similar to a Gaussian function of the quark number n at high T, and 
to a geometric series with negative n components. The Lee- Yang zero trajectory is a 
circle with two-branch at high T and two circle at low T. The eigenvalue distribution 
was almost independent of the lattice size in our simulations on 8'^ and 10'^, which 
may suggest that the canonical partition function is insensitive to the spatial volume. 
Assuming the obtained n-dependence holds for larger n, we have shown that the T- 
dependence of the canonical distribution and Lee- Yang zero trajectory are related 
to the phase transition. The canonical distribution and Lee- Yang zero trajectory 
distinguish the crossover and first order phase transition. Hence the investigation 
of the canonical partition function and the Lee- Yang zeros may provides a way to 
distinguish the phase transition and fictitious signals caused from sign problem. It 
would be interesting to combine the present approach and the ordinary techniques 
such as the finite size scaling of the Lee- Yang zero. Maxwell construction of the 
canonical formalism, etc. 

For confirmation of the results found in the present work, several points should 
be clarified further. Our analysis was performed on the small and coarse lattices 
with heavy quark mass. It is important to remove these lattice artifacts on the 
eigenvalues of the reduced matrix. The gap of the eigenvalues is sensitive to the 
quark mass. The tail of the eigenvalue distribution is sensitive to the quark mass 
and lattice spacing. These quantities are related to the pion mass and highest energy 
level, which determine the critical values of for the low temperature limits. It is 
also important to examine the A^^-scaling law for larger temporal lattice size. 

It is also important to investigate the lattice artifacts on the canonical partition 
function. In addition, the finite-size scaling and the improvement of the overlap are 
necessary task to identify the phase transition point. In particular, the large quark 
number sector of the canonical partition function plays an important role in the 
phase transition. For the improvement of the overlap, a multi-ensemble reweighting 
or histogram method may be useful.'^ 

With further confirmations, it will be interesting to study the thermodynamical 
properties of QCD at low temperatures. 

Acknow^ledgement s 

We thank S. Ejiri, Ph. de Forcrand, S. Hashimoto, M. Hanada and H. Mat- 
sufuru for valuable comments and stimulating discussions. Parts of this work were 



36 



K. Nagata, S. Motoki, Y. Nakagawa, A. Nakamura and T. Saito 



done during the stay at the workshop "New Type of Fermions on Lattice" held at 
YITP. We thank to A. Ohinishi, T. Misumi, D. Adams, C. Hoelbhng for stimulating 
discussions and valuable information. KN specially thanks to T. Hatsuda for the hos- 
pitality and encouragement during the stay at Tokyo University. AN acknowledges 
the hospitality and discussions to M. Yahiro at Kyushu university. This work was 
supported by Grants-in-Aid for Scientific Research 20340055, 20105003, 23654092 
and 20105005. The simulation was performed on NEC SX-8R at RCNP, NEC SX-9 
at CMC, Osaka University, and System A and System B at KEK. 

Appendix A 

Miscellaneous for reduction formula 

A.l. Pair nature of the eigenvalues 

We show a proof of a pair nature of the eigenvalues of the reduced matrix. The 
pair nature is a consequence of the symmetry of the reduced matrix, which was 
shown in Ref.^S Here we present a brief proof. 

Substituting Eq. ()2-6ap into the 75-hermiticity relation leads to 

r n + = (^*)^ n + (^*y')- (^-i) 

n=l n=l 

Since the 75 hermiticity holds for V/i € C, Eq. (|A-ip also holds for G C. 

Let ^0 a value of fugacity which set l.h.s of Eq. ()A-ip to be zero, i.e., ^0 = — A„. 
If l.h.s is zero, then r.h.s also must be zero. Hence, an eigenvalue for l.h.s = must 
exist for ^ = ^o- This is satisfied by A + (Co) ""^ = 0- Thus, two eigenvalues 

appears at the same time. This procedure is applied to all the eigenvalues A„, (n = 
1, • • • iVfed)- Thus, the pair nature of the eigenvalues of the transfer matrix is proved. 

A. 2. Properties of Fugacity coefficients 

Next, we consider the property of the fugacity coefficients c„ under center trans- 
formation Z3, 

[74(2;, tj) — )• a;C/4(x, tj), 3tj, Vx, (A-2) 

where ui = exp(27riA;/3), {k = ±1) is an element of Z3. In Eq. (|2-5p . Cq is invariant 
under Z3, since it does not contain temporal link variables. From Eqs. ()2-5bp and 
(j2-4cp Q transforms as Q — )• cvQ under Z3. The fermion determinant is transformed 
as 

det A{fi) = CoC^'^"''^ det{ujQ + i), 

n=0 
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A^rcd/2 

= Co (A-3) 

n=-A'rcd/2 

Then, the fugacity coefficients c„ transforms c„ — > c^w^" under the Z3 transforma- 
tions. Cn is invariant for the triality sector n = 3m(m € Z) and covariant for the 
triahty nonzero sector n = 3m + 1, 3m + 2. 



Appendix B 

Properties of canonical partition function 

B.l. RW invariance and Triality 

Consider the Roberge- Weiss transformation!^ Temporal hnks are transformed 
under Z3 

Ui{x,ti) ^ Ui{x,ti) = ujUi{x,ti), ^ti/x, (B-1) 



and the chemical potential is shifted in the imaginary direction as 

fi fj,' fi 2TTik 



(B-2) 



T T T 3 ' 

(B-3) 

which acts on the fugacity as the rotation (.^ — )• ^^oj). It is obvious from Eq. ()A-3p 
that the transformations for the link variables and chemical potential cancel. The 
fermion determinant is invariant under the RW transformation, 

det A{n, {U}) det A{n', {U}') = det A{fi, {[/}) (B-4) 

The grand partition function is also invariant under the RW transformation, 

Z{fi) ^ Z{fi') = j DU' dei A{^i' ,{UY)e-^'^^^'\ 

DU dei A{ii,{U})e~^^^^\ 

= Zi^Ji). (B-5) 

An important property on the canonical partition function is deduced from this 
invariance .f^^} 'ED Por simplicity we set /i to be pure imaginary. Expressing the 
fermion determinant as the fugacity polynomial, the grand partition function is given 
as 

» M 

Zacim) = DU CnCe^\ (B-6) 

n=-M 

and 

+^^) = [dU J2 CnCVe^^ (B-?) 

n=-M 
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Here note that the maximum quark number M is proportional to the spatial volume 
V and diverges at thermodynamical limit, which may spoils this discussion near 
phase transition points. Assuming the analiticity, Eq. (jB-5p leads to 

(l - exp(i^)) I DUcne^^ = 0. (B-B) 

Thus, we obtain 

Z{n) = j DUcne^^ = (mod(n,3) / 0). (B-9) 

Hence the canonical partition function must vanish for the triality nonzero sector. 
The above argument is applied to arbitrary number of flavors, which can be obtained 
by changing M. 

B.2. Canonical partition function for triality nonzero sectors 

We have seen the properties of the canonical partition function Eq. ()B-9p . How- 
ever, we showed that the canonical partition functions do not vanish for the triality 
nonzero sector, see Fig. [2DJ 

This small paradox is caused by the importance sampling at ^ = 0, where 
configurations for one Z3 sector are collected in the presence of the quarks. Effects 
of the non-zero triality sector was investigated in Ref.SS) 

Here, we consider this point. First, we classify the configuration space into three 
regions according to the location of the Polyakov loop (L) on the complex L plane; 

= {U\ - 7r/3 < arg(L) < 7r/3}, (B-10) 
= {U\ ^/3 < arg(L) < tt}, (B-11) 
R^ = {U\-Ti < arg(L) < -7r/3}. (B-12) 

The grand partition function is written as 

Z{ii)=([ +[ + /" ) DC/detZ\(^)e-^« (B-13) 
\Jri Jr2 JR3J 

As we have already seen, whole the Z{fj,) is invariant under the Roberge- Weiss trans- 
formation. Let us consider the transformation property of each component. Let Ij 
to be one of Z3 component of the partition function, 



Ij= DUdetA{fi)e-^^. (B-14) 
Jr, 

Under the Roberge- Weiss transformation, the Polyakov loop L is rotated in the 
complex L plane, then the configurations {U} move from one of Ri sector to another 
Therefore, Ij transforms as 

I,(/x) = /,+i(/i'),(/4 = /i)- (B-15) 
Therefore, each is not invariant under the Roberge- Weiss transformation. 
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In the importance sampling at /i = in the presence of the quarks, configurations 
for one Z3 sector are extensively collected, which cover a part of the configuration 
space, e.g. Ri. This causes the non vanishing of the triality nonzero sectors in the 
canonical partition function. 

In order to satisfy the triality is to include other Z3 sectors by using the RW 
transformation , 

Z(^)=/i(;u)+/2(/u)+/3(m) 

= + + 27rir/3) + 47riT/3) 

= j Yl 3c„re-^«, (B-16) 

n,mod(n,3)=0 

which contains only mod(n, 3) = terms. 



Appendix C 

Calculation of coefficients Cn 

The fugacity coefficients c„ of the fermion determinant can be obtained in the 
reduction formula, 

n = -Afred/2 

Their values vary from order one to order 10^'^'' even on the small 4^ lattice. They 
cannot be handled in the double precision. 

The simplest way is to use an arbitrary accuracy libraries. This is more than 
necessary. To express Ck, we need wide range of floating numbers, but we do not 
need very high precision. In other words, we need wide range of the exponent, but 
we do not need very huge significant numbers. 

We express each real and imaginary parts of in a form of 

a X 6^, (C-2) 

where 

1 < |a| < b, (C-3) 

and a is a double precision real and L is an integer. We employ the "module" in 
Fortran 90, which allows us to define a new type of data and mathematical operation 
among them. The base b can be any number, and we set it to be 8. 

There are several way to get in Eq. (|CTp . The simplest way is to use in a 
recursive way: 

M M-1 

Y,c'ke = {Bo+Bii)Y,Cke (c-4) 

fc=0 fc=0 
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and 

Cq = BqCo 

C'k = Bk-iCk + BkCk-i (A; = 1, 2, • • • , M - 1) 
C'm = BiCm-^1 (C-5) 

We calculate several cases by this method, and by a high accuracy library, FMLIB.!^ 
We got the same results. 

A smarter way is a divide-and-conquer method. See ([1]). Here for simplicity we 
assume = 2^^, but this can be loosened c cl x c2 is an operation to determine 
c from ci and C2 where 

(c(0) + c(l) * X + ... + c(iV3) * x^) ^ 

(ci(0) + ci(l) * X + ... + ci{Ni) * x^i) (C-6) 

X(C2(0) + C2(l) * X + ... + C2(iV2) * X^2) (C-7) 

and iVg = iVi + iV2. 



Algorithm 1 Divide-Conquer calculation for the coefficients c„ for nj=i('^(^) + ^) ~ 

Input a, Output c) 
Set N 
read a 

CALL DandQ ( a, c, N) 

Recursive SUBROUTINE DandQ ( a, c, N) 

IF N corresponds the Bottom, RETURN 

CALL DandQ (a(l), cl, N/2) 

CALL DandQ (a(N/2+l), c2, N/2) 

c ^ cl X c2 

RETURN 

END SUBROUTINE 



When is large (i.e., we are near to the "top" of the recursive level) , Eq. ()C-7p 
is easy to vectorize or parallelize. For small A^ (i.e., we are near to the "bottom" of 
the recursive level) , we handle many calculations of the type of Eq. ()C-7p , then it is 
also easy to vectorize or parallelize. 
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